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NORMAL  AND  RADIAL  IMPACT  OF  COMPOSITES 
WITH  EMBEDDED  PENNY-SHAPED  CRACKS 

by 

G.  C.  Sih 

Institute  of  Fracture  and  Solid  Mechanics 
Lehigh  University 
Bethlehem,  Pennsylvania  18015 

and 

* 

E.  P.  Chen 
Sandia  Laboratories 
Albuquerque,  New  Mexico  87115 

ABSTRACT 

A  method  is  developed  for  the  dynamic  stress  analysis  of  a  layered  composite 
containing  an  embedded  penny-shaped  crack  and  subjected  to  normal  and  radial  im¬ 
pact.  The  material  properties  of  the  layers  are  chosen  such  that  the  crack  lies 
in  a  layer  of  matrix  material  while  the  surrounding  material  possesses  the  aver¬ 
age  elastic  properties  of  a  two-phase  medium  consisting  of  a  large  number  of 
fibers  embedded  in  the  matrix.  Quantitatively,  the  time-dependent  stresses  near 
the  crack  border  can  be  described  by  the  dynamic  stress  intensity  factors.  Their 
magnitude  depends  on  time,  on  the  material  properties  of  the  composite  and  on 
the  relative  size  of  the  crack  compared  to  the  composite  local  geometry.  Re¬ 
sults  obtained  show  that,  for  the  same  material  properties  and  geometry  of  the 
composite,  the  dynamic  stress  intensity  factors  for  an  embedded  (penny-shaped) 
crack  reach  their  peak  values  within  a  shorter  period  of  time  and  with  a  lower 
magnitude  than  the  corresponding  dynamic  stress  factors  for  a  through-crack. 


*This  work  was  completed  when  Dr.  Chen  was  a  faculty  member  at  Lehigh  University. 
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INTRODUCTION 


Advanced  composite  materials  are  multi -phased  nonhomogeneous  materials  with 
anisotropic  properties.  This  complicates  the  stress  analysis  for  fracture,  par¬ 
ticularly  if  the  loading  is  time-dependent  and  the  geometry  involves  sharp  edges 
such  as  a  crack.  As  a  result,  conventional  and  mathematical  techniques  for  dy¬ 
namic  fracture  generally  fail  to  yield  accurate  results. 

An  effective  approach  for  finding  dynamic  stresses  in  a  nonhomogeneous  com¬ 
posite  containing  a  through  crack  has  been  developed  [1]  by  utilizing  both  the 
Laplace  and  Fourier  transforms.  The  transient  boundary,  symmetry  and  continuity 
conditions  were  formulated  by  integral  representations  in  terms  of  the  rectangu¬ 
lar  Cartesian  coordinates  x  and  y  and  the  results  for  the  stress  intensity  fac¬ 
tors  are  determined  numerically  by  solving  a  standard  integral  equation  in  the 
Laplace  transform  plane.  The  crack  geometry  was  assumed  to  be  extended  infinitely 
in  the  z-direction  or  through  the  side  wall  of  the  composite  specimen.  Many  of 
the  failures  in  composites,  hov/ever,  were  observed  [2]  to  initiate  from  embedded 
mechanical  imperfections  such  as  air  bubbles,  voids  or  cavities.  Hence,  a  more 
realistic  modeling  of  the  actual  flaw  geometry  would  be  an  embedded  crack  that 
has  finite  dimensions  in  all  directions.  This  immediately  suggests  a  three-di¬ 
mensional  elastodynamic  crack  problem  which  cannot  be  solved  effectively  by  ana¬ 
lytical  means  unless  symmetry  prevails.  One  approach  for  obtaining  a  solution  is 
to  extend  the  integral  transform  formulation  for  a  through  crack  in  rectangular 
coordinates  [1]  to  that  of  an  embedded  crack  in  cylindrical  polar  coordinates. 

This  necessitates  the  use  of  Hankel  transforms  instead  of  Fourier  transforms. 

Although  no  attempt  will  be  made  to  analyze  the  failure  of  the  composite  due 
to  impact,  the  dynamic  stress  intensity  factors  k^(t)  and  k2(t)  can  be  readily 
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used  in  a  given  fracture  criterion,  say  the  strain  energy  density  theory  [3], 
for  determining  the  allowable  level  of  impact  load.  The  new  results  can  also 
assist  the  construction  of  composite  materials  for  establishing  impact  tolerance. 
In  this  case,  failure  is  assumed  to  initiate  from  a  damage  zone  of  material  in 
the  composite  that  can  be  approximated  by  an  embedded  crack.  The  time-dependent 
characteristics  of  the  stresses  for  the  through  and  embedded  crack  geometries  are 
compared  and  studied  for  different  elastic  properties  and  dimensions  of  the  com¬ 
posite.  In  particular,  the  phenomenon  of  elastic  waves  reflecting  from  the  crack 
to  the  interfaces  within  the  composite  can  be  exhibited  numerically  when  their 
neighboring  boundaries  are  sufficiently  close  to  one  another.  As  time  becomes 
very  large,  all  of  the  results  in  this  report  reduce  to  the  corresponding  static 
solutions  [4]. 

AXIAL  SYMMETRIC  DEFORMATION:  PENNY-SHAPED  CRACK 

Consider  a  penny-shaped  crack  of  radius  a  that  lies  in  a  layer  of  material 
of  thickness  2b  with  material  properties  ,  v-j ,  p-j.  This  layer  is  bonded  be¬ 
tween  two  media  with  properties  ^2’  ^2  illustrated  in  Figure  1.  With 
reference  to  the  system  of  coordinates  (x,y,z),  the  z-axis  coincides  with  the 
center  of  the  crack  and  is  normal  to  the  crack  situated  in  the  xy-plane.  The 
outer  boundaries  of  the  composite  are  assumed  to  be  sufficiently  far  away  from 
the  crack  such  that  the  reflected  waves  will  have  a  negligible  influence  on  the 
local  stresses.  Only  those  impact  loads  that  produce  an  axisymmetric  wave  pat¬ 
tern  will  be  considered. 

For  an  axially  syiranetric  deformation  field,  material  elements  are  displaced 
only  in  the  radial  and  axial  direction  and  remain  unchanged  in  the  9-direction. 
With  reference  to  the  cylindrical  polar  coordinates  (r,9,z)  in  Figure  1,  the 
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two  nonzero  displacement  components  can  be  expressed  in  terms  of  the  wave  poten¬ 


tials  (f).(r,z,t)  and  Tp.(r,z,t)  as  follows; 

J  J 


3(J). 

f  u  )  —  — ^  ^ 

^^r^.  3r  3z 

\J 


3({)  ij;. 

( 11  )  =  - i  i  ^ 

^“z^.  3z  3r  r 

vJ 


where  j  =  1  refers  to  the  layer  with  the  crack  and  j  =  2  to  the  surrounding  ma¬ 
terial.  The  four  nontrivial  stress  components  are  given  by 


3(J).  3iiJ. 

^ ^ _ J,  _  _ iA  +  A 


(a^)  -  9^  (gp  -  92^ 

J 


,  3({).  3i|j. 

(Oj)  =  2uj  7  -  35-)  + 


94) .  dip  •  ip  • 

(Oz)  _  -  2yj  +  — )  +  XjV  (})j 

.  9c|).  dilJ.  .  3(J). 

(Vz>.  =  I'j  C|i-  (2  ^  ^ 


in  which  A-  and  y.  are  the  Lam^  constants  and  represents  the  operator 
J  J 


3^  +  i  1_  +  ll 

3r^  r  3r  3z‘ 


The  governing  equations  can  thus  be  obtained  from  the  equations  of  motion  which 


yield 


(3) 


^  r  3r  ^  3z^  c|T  3t^ 

3F^‘^?3F“ 


with  c-jj  and  C2j  being  the  dilatational  and  shear  wave  speeds: 

'ij  =  ’  '=20  “ 

If  the  composite  body  is  initially  at  rest,  the  Laplace  transform  of  equations 
(3)  further  give 

<\2x*  ‘ij.* 


iii  +  i  !!i  +  i!l .  hL  ** 

3r^  *  f  oT  *  ^  C?7  ♦  j 


★  ★ 


(,  .  aip.  ij).  3^t|>.  2  * 


Here,  p  is  the  transform  variable  in  the  Laplace  transform  pair: 


f  (p)  =  /  fit)  exp(-pt)dt 
0 


fit)  =  i  /  f*Cp)  exp(pt)dp 

ZTTl 


The  abbreviation  Br  stands  for  the  Bromwich  path  of  integration.  Moreover,  since 
the  composite  geometry  is  symmetrical  about  the  xy-plane,  it  suffices  to  consider 
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only  the  solution  in  the  upper  half-space,  z  ^  0.  For  the  penny-shape  crack  ge¬ 
ometry,  the  Hankel  transform  pair  [5]  may  be  used: 


f*^(s)  =  /  xf(x)  J^(sx)clx 
0 

f(x)  =  /  sf'’ts)  J„Csx)ds 
0 


(7) 


where  J  is  the  nth  order  Bessel  function  of  the  first  kind.  Applying  equations 
n 

(7)  to  (5),  the  following  results  are  obtained: 


(j)^(r,z,p)  =  /  [A^^^(s,p)e  +  A^^^(s,p)e  ]  JQ(rs)ds 

'  0 

(8) 

il^tCr.z,p)  =  /  [B^^Hs,p)e  +  B^^Hs,p)e^^^  ]  J^(rs)ds 
0 


for  the  cracked  layer  and 


4>2('">2»p)  *  /  C^^^s,p)e  JQ(rs)ds 

0 

i|;*{r,z,p)  =  /  C^^^.s,p)e  J^(rs)ds 

0 

for  the  surrounding  material .  The  quantities  y^.^.  are  given  by 

2  1/2  „2  1/2 

^1  j  "  ’  ^2j  =  ^ 


C9) 
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The  six  unknowns  .  are  determined  from  a  given  set  of  tran- 

sisnt  boundary j  synmetry  and  continuity  conditions. 

NORMAL  IMPACT 

Let  the  penny-shaped  crack  be  subjected  to  a  uniform  impact  load  such  that 
the  upper  and  lower  surface  will  move  in  the  opposite  direction.  The  magnitude 
of  this  normal  load  is  and  since  it  is  applied  suddenly  from  t  *  0  and  main¬ 
tained  at  a  constant  value  thereafter,  the  Heaviside  unit  step  function,  H(t), 
will  be  used,  i.e.,  -aQH(t).  Making  use  of  equations  (6),  the  conditions  on  the 
plane  z  =  0  for  r  _<  a  and  r  ^  a  take  the  forms 


(a*)  tr,o,p)  =  -  CT*3)^Cr,o,p)  =  0,  0<r<a 


(u!)  (r,o,p)  =  0;  (t*_)  (r,o,p)  =  0,  r  >  a 


If  the  interfaces  at  z  =  ±b  is  bonded  perfectly,  the  stresses  and  displacements 
can  then  be  considered  continuous  across  these  planes,  i.e., 

(a*)  {r,b,p)  =  (a  )  {r,b,p) 

2  1  ^2 

(12) 

(t*  )  (r,b,p)  =  (x*  )  (r,b,p) 


*Thpre  is  no  loss  in  generality  in  formulating  the  problem  in  terms  of  a  unifoim 
step  load.  The  principle  of  superposition  may  be  used  to  obtain  the  solution  for 
general  loading  from  a  series  of  step  loading  solutions  as  discussed  in  [IJ. 
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and 


(u  )  (r,b,p) 
1 


(u^)  Cr,b,p) 
2 


'k  k 

(u  )  (r,b,p)  =  (u  )  (r,b,p) 
z  1  z  2 


(13) 


Under  these  considerations,  the  six  functions  may  be  ex¬ 

pressed  in  terms  of  a  single  unknown  A(s,p)  as  indicated  by  equations  (A.l)  in 
the  Appendix. 


F^exLkotm  Znte.gACLl  zqmtloYU>.  Without  going  into  details,  the  function  A{s,p) 
can  be  obtained  from  the  system  of  dual  integral  equations 


00 

/  A(s,p)  jQ(rs)ds  =  0,  r  >  a 
0 


00  a 

J  5Pj{s,p)  A(s.p)  J„Crs)ds  =  -  r  <  a 


(14) 


in  which  Pj(s,p)  is  a  known  function: 

2  "2(y 

^  sAj(1-kP  ^^^^21^  ■  ^^^11^21^*^*^^^^  ~ 


■'^21 


)b 


] 


(15) 
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The  form  of  A(s,p)  that  satisfies  equations  (14)  can  be  found  from  Copson  [6]: 


—  a  a^/^  1 

2v,ptl“<fT  I  J,/2(5a5)<i5 


yp 

Here,  Jy2  Bessel  function  of  the  first  kind  and  Aj(?,p)  satis¬ 

fies  the  Fredholm  integral  equation 


At(?>P)  +  I  AT(ri>p)  M,(5,n,p)dn  =  5 
^  0  ^  ^ 


whose  kernel 


MjCc.n.p)  =  /Cn  /  s[Pj(f,p)  -  1] 


1/  tPT(f’P)  ■  sin(sS)  sin(sn)ds 

It  1  la 
0 


is  symmetric  in  5  and  n.  Figures  2  to  4  show  the  numerical  results  of  equation 

(17)  by  varying  ii2/Pi  while  p-j  =  P2  and  v-j  =  V2  =  0.29  are  kept  the  same 

for  all  cases.  The  function  A*(5,p)  evaluated  at  the  crack  border,  C  =  1.  gov- 

* 

erns  the  contribution  of  the  geometric  and  material  parameters  on  ki(p)  which 
represents  the  Laplace  transform  of  the  stress  intensity  factor. 

*  .  . 

Znt&n&'ity  IclcX-ok  ^oa.  noAmaZ  impacit.  In  order  to  evaluate  k-jCp)  or 
k^(t),  the  stresses  in  the  matrix  layer  are  first  expanded  in  terms  of  the  local 
coordinates  r^  and  0.j  for  small  values  of  r-j.  The  local  coordinates  (r^,e^)  are 
related  to  (r,e)  in  Figure  1  as  follows: 


a  +  cos0^  =  r  cose 


(19) 


sin0^  =  r  sine 


The  leading  term  in  the  Laplace  transform  of  the  local  stresses  that  possess  the 
IZ/rj"  singularity  is 


k-]  (p) 


p  ir  0 


(20) 


Application  of  the  Laplace  inversion  theorem  yields  the  dynamic  stress  field 
around  the  crack  border  as  a  function  of  time.  The  result  is 


k,(t)  0,  0n  30, 

(^v’)  (r,,0,  .t)  =  — - cos  ^  (1  -  sin  sin  -j-)  +  0(rp 

rill  ^  d  d  d  i 

ki(t)  0,  ^  ^ 

(oq)^  (r^,ei,t)  =  - - 2v.|  cos  ^  +  0(rp 


0,  30, 


k-i  (t)  0,  U-I  -J^l 

{o)  (r,  ,0,  ,t)  =  — - cos  yL  (1  +  sin  j-  sin  -«-)  +  0(rp 

z  1  1  1  ^  d  d  d  \ 

k,(t)  0,  0,  30, 

(t^z)  (>^1  ’^i^t)  =  — - cos  sin  ^  cos  —  +  0(rp 

1  /2r.| 


(21) 


and  k-|(t)  becomes 


2a  /a  ,  At(I.p)  q. 


Br 


(22) 


Note  that  equation  (20)  is,  in  fact,  the  Laplace  transform  of  equation  (22). 
Hence,  the  functional  dependence  of  r.|  and  0i  is  not  affected  by  the  Laplace 
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transformation  and  can  be  evaluated  separately.  This  observation  was  first  made 
by  Sih,  Ravera  and  Embley  [7]. 

'it 

Making  use  of  the  results  for  Aj(l,p)  in  Figures  2  to  4,  k.|(t)  in  equation 
(22)  can  be  found  as  given  in  Figures  5  to  7.  The  dynamic  stress  intensity 
factors  k.|(t)  for  the  penny-shaped  crack  exhibit  an  oscillatory  behavior  rising 
quickly  to  a  peak.  As  time  increases,  all  curves  approach  the  static  value  of 
k-j  =  2a^»^/Tr  [4].  For  a  crack  diameter  to  layer  thickness  ratio  of  a/b  =  1,  the 
peaks  of  the  k-j(t)  curve  are  sensitive  to  changes  in  the  shear  moduli  ratio 
Figure  5  indicates  that  k.,(t)  tends  to  decrease  in  amplitude  as  is  reduced 

from  0.1  to  10.0.  The  influence  of  the  composite  interface  on  k.|(t)  is  exhibited 
in  Figures  6  to  7.  When  the  shear  modulus  of  the  surrounding  material  U2  is 
much  smaller  than  the  matrix  layer  with  the  dynamic  crack  border  stress  in¬ 
tensity  increases  as  the  crack  diameter  becomes  large  in  comparison  with  the  layer 
thickness.  This  effect  is  clearly  evidenced  in  Figure  6.  As  expected,  k-^{t)  in¬ 
creases  with  decreasing  a/b  when  the  shear  modulus  of  the  cracked  layer  is  made 
smaller  than  the  surrounding  material,  i.e.,  <  y2  as  illustrated  in  Figure  7. 

The  result  of  Embley  and  Sih  [8]  is  recovered  for  the  homogeneous  case,  y.|  =  y2. 

RADIAL  IMPACT 

If  the  penny-shaped  crack  is  sheared  uniformly  in  the  radial  direction  such 
that  axial  symmetry  is  preserved,  then  (J>j. 0)^,2, p)  and  ii^j(r,z,p)  in  equations  (8) 
and  C9)  remain  valid.  Let  this  shear  of  magnitude  be  applied  suddenly  and 
hence  the  surface  tractions,  -x^HCt),  are  to  be  specified  for  0<r<a  with  H(t) 
being  the  Heaviside  unit  step  function.  Laplace  transform  of  the  conditions  on 
the  plane  z  =  0  thus  become 


-11- 


(4)^(r.o.p)  =  -^ 


TT;  (<^y)  Cr>o,p)  =  0,  0<r<a 
P  z,  ^ 


T^P 

(Uy,)  (r,o,p)  =  0;  (a  )  Cr,o,p)  =  0,  r  >  a 
r  1  z  ^ 


Continuity  of  the  stresses  across  the  interface  z  =  b  is  satisfied  if 


(a  )  (r,b,p)  =  (a  )  (r,b,p) 
z  1  z  2 


(<y*)  (r,b,p)  =  (a*  )  (r,b,p) 


and  the  same  requirement  is  imposed  on  the  displacements: 


it  if 

(Uy,)  (i^.b.p)  =  (u^)  {r,b,p) 


if  if 

(u,)  Cr,b,p)  =  Cu,)  (r,b,p) 
z  1  z  2 


JntzgA.aJL  zquoZioru,.  As  in  the  case  of  normal  impact,  the  six  unknown  func¬ 
tions  A^^^(s,p),  A^^^ (s,p) , . . . ,  C^^^(s,p)  in  equations  (8)  and  (9)  can  be  ex¬ 
pressed  in  terms  of  a  single  unknown  B(s,p).  Refer  to  equations  (A. 5)  in  the 
Appendix.  Hence,  equations  (24)  and  (25)  are  satisfied.  The  remaining  boundary 
conditions  in  equations  (23)  are  employed  to  obtain  the  system  of  dual  integral 
equations 


/  B(s,p)  J.J (rs)ds  =  0,  r  >  a 


/  sPjj(s,p)  B(s,p)  J^(rs)ds 


r  <  a 


in  which 

Pll(s,p)  = PjCs.p)  (27) 

where  Pj(s,p)  is  already  known  through  equation  (15)  while  Aj(s,p)  and  Ajj(s,p) 
are  given  by  equations  (A. 2)  and  (A. 6),  respectively. 

Solving  for  B(s,p)  [6],  it  can  be  shown  that 


B(s,p)  =  - 


2  4v^pC1-k|)  J 


/  yf  AttC^.P)  Jv?(saC)d5 


and  AttC^.p)  satisfies  the  Fredholm  integral  equation  of  the  second  kind: 


Att(C.p)  +  j  Ar.Cn.p)  M,T(?,n.p)dTi  =  Z 


whose  kernel  takes  the  form 


MTT(?.n,p)  =  ^  /  s[Ptt(|-,  p)  -  1]  J-5/o(s5)  J3/2(sn)ds 


Plots  of  Ajj(l,p)  as  a  function  of  C2i/pa  are  shown  in  Figures  8  to  10  for  dif- 
ferent  values  of  v^2^^'\  a/h.  The  curves  show  that  AjjCl,p)  rises  rapidly  at 


first  and  then  levels  off. 
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Stn.2M&  iacXofi  ion.  fiadial  mpaxit.  The  dynamic  crack  border  stress 

field  corresponding  to  radial  shear  can  be  obtained  in  the  same  way  and  expressed 
in  terms  of  the  coordinates  (r^,0^)  in  equations  (19): 


k.(t)  6,  0,  30, 

(o'y,)  (r,  ,0T,t)  =  - sin  (2  +  cos  cos  -j^)  +  0(r?) 

1  '  '  yzr^  ^  ^  1 


k2(t)  0, 

(r,.0,,t)  =  ^ - 2v,  sin2^+0(rp 

Bill  ' 


kpCt)  0,  0,  30, 

((J,)  (r,  ,0,  ,t)  - - sin  cos  ^  cos  -5-  +  0(r?) 

z  1  I  1  ^  z  z  z  I 


kpCt)  0,  0,  30, 

(Xrz)  C^i.e-i.t)  =  - - cos  ^  (1  -  sin  j-  sin  -^)  +  0(rp 

1  y  2r  1 


Note  that  k2(t)  can  be  evaluated  from 


I  rti  -  I  ptj„ 


once  Ajj(l,p)  as  given  by  Figures  8  to  10  is  known. 

The  numerical  results  in  Figures  11  to  13  for  k2(t)  as  a  function  of  time 
refer  to  =  p2  and  v-j  =  V2  =  0.29.  The  curve  with  li-j  =  ^2  ''s  the  solution  for 
the  homogeneous  material  treated  previously  by  Embley  and  Sih  [8].  In  general, 
k2{t)  oscillates  with  time  and  can  be  greater  or  smaller  than  the  corresponding 
homogeneous  solution  depending  on  whether  iJ2/l->-|  ^  O'"  ^  Figure  11 

displays  the  variations  of  k2Ct)  for  different  values  of  ^2/^!]  whii®  s/t  is 
fixed  at  unity.  The  influence  of  the  ratio  of  crack  size  with  layer  thickness 


is  exhibited  in  Figures  12  and  13  for  Vi2^^1  ~  lO-O*  respectively. 

These  two  cases  show  the  opposite  effect  which  is  to  be  expected. 

CONCLUDING  REMARKS 

The  previous  discussion  has  shown  that  the  dynamic  stress  intensity  factors 
for  an  embedded  crack  can  be  evaluated  analytically  by  a  method  similar  to  that 
developed  for  a  through  crack  [1].  An  important  consideration  is  to  compare  the 
results  for  these  two  crack  configurations  and  to  draw  some  general  conclusions. 
First  of  all,  the  k^(t)  or  k2(.t)  factor  for  the  penny-shaped  crack  tends  to  rise 
more  quickly  than  the  through  crack,  i.e.,  the  peak  value  of  k^(t)  or  k2(t)  is 
reached  within  a  shorter  period  of  time.  This  is  because  waves  emanating  from 
the  neighboring  points  on  the  periphery  of  the  penny-shaped  crack  interfere  with 
each  other  much  earlier  as  compared  to  a  line  (or  plane)  crack  where  the  waves 
must  travel  from  one  end  to  the  other  before  interference  can  take  place.  In 
general,  the  maximum  value  of  k-j(t)  or  k2Ct)  for  an  embedded  crack  is  lower  than 
that  for  a  through  crack.  For  example.  Figure  5  gives  a  peak  value  of  approxi¬ 
mately  1.6  for  TTk-j  (t)/2aQj/a  which  corresponds  to  a/b  =  1.0  and  =  O-l-  This 

occurs  at  C2it/a  1.6  and  yields  k^(.t)  -  1.02  The  corresponding  case  of  a 

through  crack  [1]  renders  k^(t)  -  2.40  o’qV^  3nd  C2it/a  =  3.0.  The  difference  in 
k-j  (t)  is  more  than  a  factor  of  two  and  is  more  pronounced  as  the  ratio  a/b  is  in¬ 
creased.  For  embedded  cracks  that  are  non-circular  in  shape,  approximate  esti¬ 
mates  of  k^(t)  can  be  made  by  taking  the  solution  for  the  through  crack  as  an 
upper  limit  and  that  of  the  circular  crack  as  a  lower  limit. 

In  the  absence  of  axi symmetry,  the  dynamic  stress  analysis  will  become  ex¬ 
ceedingly  difficult  and  it  will  be  more  feasible  to  solve  the  crack  problem  nu¬ 
merically.  In  such  cases,  the  solutions  obtained  here  can  be  used  to  guide  the 
development  of  numerical  procedures. 
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APPENDIX:  EXPRESSIONS  FOR  ^  (s,p) , . . . ,  C^^^s.p) 

UoHxnaZ -impact.  The  functions  A^'*^(s,p),  A^^^  (s  ,p) , . . . ,  C^^^(s,p)  for  the 
wave  potentials  in  equations  (8)  and  (9)  can  be  expressed  in  terms  of  a  single 
unknown  A(s,p)  for  normal  impact 


A<'>(s,p)  .  [1 


6W)e' 


■^'^21'’, 


SYiiO 


■(y„+y 


21 


)b 


A(s,p) 

"l 


A<^'Cs,p)  =  -  [sY„e  (s“+Y2i)e'^^"'' 


ACspP) 


b(')(3,p)  _  ,(2)j(2)Jll'>j 


b(2)(s.p)  =- 

-  s(Y2TY22)B^^^e  +  s (Y2i+Y22)B^^^e^^^ 

'  s^^i2Y22  '^®<d2-Yn>''‘'’®  ■"  =(Yn+Y,2)e^" 

Ms^-Y2,Y,2)B('>e'"2''Ms^Y,2,Y,2)B<2)e"21^ 
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in  which  Aj  stands  for 


A.(s.p)=4y„ 


+  6 


(l)e-'"ll^ 


and  5^^^  6^^\...,  are  further  expressed  1n  terms  of  e^'\  e^^\..., 
as  the  following: 

6(')(s,p)  =  (e<’)e<«>  -  e'2>e(^>)/Ce(’>e'«>  -  e^^le'^)) 

,(2),s,p)  =  (e(«)e(«>  -  e(2)e(8))/(e(”e<«>  - 

d<3)(s,p)  >  (e<l>e(^>  -  e‘3)e(5))/(a(l)e(«)  -  e<2)e(3)) 

d(^)(s,p)=(e<'>e(8)-e(^>e'3))/(e<'>e<«>-e(3)e(5)) 


(A.2) 

e(8) 


(A.3) 


The  quantities  in  equations  (A. 3)  are  complicated  functions  of  the  materials  pa¬ 
rameters  and  transform  variables.  They  are  given  by 

e^^hs,p)  =  -  sy2^  +  y^(s"-Y^2Y22')'  ^22(5^-^21^12)^ 

e(2)(s.p)  =  SY2T  -  p^(s^-yi2T2?  (^21+^22) (^'^22)  ’  ^22 (^'''^21^1 2)^ 

e^^)(s,p)  =  I  (s2+y|i)  -  y^(s^-Y^2Y22)  (^^'^^22) ^^"'^1 1^22) 

+  s^22(YirYi2)] 
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=  1  (S^^+Y^,)  -  y^(s^-Yi2Y2P'  (s"-^22)(s"-^11Y22) 


s^Y22(Yn+Ti2)] 


e<®’(s,p)  -  -  1  (s^n|,)  +  p^  (s^-Yi'gYj;)'  '^='^12'^2r^22> 


(A. 4) 


e^^)(s,p)  =  -  1  Cs^+yfi) 


'nz'fzz' 


i2(Y2i+T22) 


^  (s^+Y22)(s^+Y2iY-i2)^ 


e^^^s.p)  -  ii^Cs^-y^2Y22^  '^^12^^^"^n‘^22^  ^  2  ^^^'^^22^ 


e^^^(s,p)  =  -  sy 


1 1  U*!  (^^"^*12^22 ^ 


Cyi2(s^+yny22)  ■  f  (s^+y22)(Yii+y-]2)] 


RacUaZ 'impact.  For  radial  impact,  A^^^(s,p),  A^^^  (s,p) , . . . ,  C^^^(s,p)  in 
equations  (8)  and  C9)  can  be  expressed  in  terms  of  B(s,p)  as 


A<')(P,P)  =  -  [sy2,(««>  -  Y  I 


(A. 5) 


a(2)rc  Tc  a:(3)  “  ^^21*^^  .1  /  2 .  2  x  "^^ll'^21 

A^  ^(s,p)  =  [sy5,e  (6^  '  -  6^  'e  )  +  y  (s^+yoi)e  ] 
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where 


‘II 


=  _E 


Yj,  [5 


C2) 


+ 


.  ,(4)  -2''2i'> 


-  5 


(1) 


b 

] 


(A. 6) 


The  remaining  functions  B^^^(s,p),  B^^^Cs.p),  etc.,  can  be  related  to  B(s,p) 
through  A^^Us.p)  and  A^^^(s,p)  since  the  last  four  expressions  in  equations 
(A.l)  for  normal  impact  also  apply  to  radial  impact. 
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Normal 
Impact  o-(t) 


Figure  1  -  Penny-shaped  crack  embedded  in  a  matrix  layer 
under  normal  and  radial  impact 
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7rk|  (t)/2cro*/a” 


Figure  5  -  Dynamic  stress  intensity  factor  k-|(t)  for  penny-shaped 
crack  with  a/b  =1.0 
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Figure  6  -  Dynamic  stress  intensity  factor  k^(t)  for  penny-shaped 
crack  with  y2/^]  =0.1 
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TTk,  (t) 


0  2.0  4.0  6.0  8.0  10.0 


C  2 1  t/a 

Figure  7  -  Dynamic  stress  intensity  factor  k, (t)  for  penny-shaped 
crack  with  ^2/^1  ~ 
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a/b  =  2.0 


1.0  2.0  3.0  4.0  5.0 


^21^  PQ 

Figure  9  -  Variations  of  Ajj.(l,p)  with  Cgi/pa  for.y2/i^l 
varying  a/b 


(d’lrv 


Figure  10  -  Variations  of  Ajj(IjP)  with  C2i/pa  for  -  10  and  varying  a/b 
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Figure  11  -  Stress  intensity  factor  k2(t)  versus  time  for  a  penny-shaped 
crack  with  a/b  =  1.0 
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COMPUTER  PROGRAMS 


k)djxl  i^pa-ct 


PROGRAM  seta  (INPUT»OUTPUTtPUNCH, PLOT *TAPE  99  =  PL0T) 
REAL  NON(A) »F(4»4»1) »G(4»4) tO(4) »PT (4» 

REAL  B(4)»C(4) 

REAL  LPCSO)  ,0TA(5g) 

EQUIVALENCE  (NON,B) 

COMMON  K1,K2,K3,K4  ^  ^ 

COMMON/ AliX /H  ,PtPKltPK2,  BMU,  X,  Y 

LP(1)=0.0 

OTA(1)  =  0.0 

READ  2,K1»K2»K3,K4 

K1  =°OROER^CF  SYSTEM  OF  EQUATIONS 
K2  =  NO.  OF  DISTINCT  KERNELS 
K3  =  NO,  OF  DATA  POINTS 
K4  =  NO.  OF  DATA  SETS  TO  BE  EVALUATED 
SET  UP  DATA  POINTS 
AK=K3 

DO  5  N=1»K3 
AN  =  N 

5  PT{N)=AN/AK, 

SET  UP  INTEGRATION  MATRIX 
M=K3-2 
N=K3-1 
A=K3 

A=1./(3.'*A) 

DO  10  K=2,M,2 
10  D{K)=2.*A 

DO  15  K=1,N,2 
15  D(K1=4.*A 

D(K3) =A 

CALCULATE  NONHOMOGENEOUS  TERMS 
RHS=1.0 
DO  22  1=1, K2 
PRINT  9 
9  FORMAT(lHl) 

READ  61,eMU 
61  FORMAKFIO.  5) 

DO  999  11=1, KA 
CO  35  N=1,K3 
35  NON(N)=RHS*PT(N> 
calculate  kernel  MATRICES 
CALL  CONST(I) 

DO  20  N=1,K3 
DO  20  M=1,K3 
IF<M-N) 25,30,30 
25  FCM ,N,I)=F (N,M, I) 

GO  TO  20 

30  F(M,N,I)=FU(I,PT (M) ,FT(N)) 

20  CONTINUE 

CALL  CHANGE CF,G,D, I) 

CALL  LIN£Q(G,B,C,  K3) 

DO  40  L=1,K3  ,  ^ 

PRINT  6,FTCL),N0NCL> 

6  FORHATt5X,F8.4,F15.6) 

40  CONTINUE 

LP{II+l)=NONCK3) 

OTA  (II+1)=P 

punches!,  (DTA(IX)  ,LP  (IX)  ,IX=1,19) 

CALL*LAPINV(OTA,LP) 

22  CONTINUE 
END 


FUNCTION  SIMF(I,A,8)  ^  ^ 

COMMON/ AUX /H,P,PK 1,PK 2, BMU,X,Y 
OEL=0.25^(B-A) 

IF (DEL)4G,45,50 
45  SIMP=0.0 
RETURN 

50  CONTINUE 

SA=Z(I,  A)+2  (I,B) 

S3=Z(I, A+2.*DEL) 
SC=Z{I,A+0EL)+Z(I,A+3.*DEL) 

-34- 


Sl=(0EL/3. )»(SA+2.»SB+4.*SC) 

IF(S1.£C.O.O)  GO  TO  45 
K=8 

55  S3  =  SB+SC 

D£L=0.5»0EL 
SC=Z(I,  A  +  DEL) 

J=K-1 

00  5  N=3,J»2 
AN=N 

)  SC=SC«-Z«I,AfAN*OEL) 

S2=(0EL/3.>*CSA*2.*SB+4.*SC) 

DIF=AaSf  (S2-S1) /S13 
ER=0.01 

IF(DIF-£R>  30t25,25 
JO  SIMP=S2 
RETURN 
>5  K=2»K 

S1  =  S2 

IF  (K-2048)  35,35t40 

*42  FORHATtlxJi^iNT.  DOES  NOT  CONVERGE  *,I3,2F9.4) 
PRINT  60,X,Y 
jO  FORMAT(2F10.5) 

00  70  J=l,10 
DIP=J 

OIP=OIP/10. 

H=Za,OIF) 

PRINT  60, W 
3  CONTINUE 
CALL  EXIT 
ENO 


7 

7 

7 

10 

11 

24 

30 
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40 

41 


SUBROUTINE  CHANGE (F,G, 0,1) 
REAL  F(4,4, 1) ,Gf4,4)  ,0(4) 
COMMON  K1,K2,K3,K4 
00  10  N=1,K3 
00  10  M=l,<3 
G(M,N)  =F(M,N,I) *0(N) 

10  CONTINUE 

00  20  N=1,K3 
20  G(N,N)=G(N,N)+1.0 

RETURN 
ENO 


SUBROUTINE  LINEGCA ,8,T,N) 
REAL  A(N,N)  ,E(N)  ,T(N) 

DO  5  1=2, N 

5  A(I,1)  =  A(1, 1)/A(1, 1) 

00  10  K=2,N 
M=K-1 

00  15  1=1, N 
15  T(I)=A(I,K) 

00  20  J=lfM 
A(  J,K)  =TtJ) 

J1=J+1 

00  20  I=J1,N 
T(I)=T(I)-A  CI,J)*A{J,K) 

20  CONTINUE 

A(K,K)=T{K) 

IF(K.EQ.N)  GO  TO  10 
M=K  +  1 

00  25  I=M,N 
25  A(I,K)=T(I)/A(K,K) 

10  CONTINUE 
*  BACK  SUBSTITUTE 
00  30  1=1, N 
T(I)=8(I) 

M=I  +  1 

IFtM.GT.N)  GO  TO  30 

00  30  J=M,N 

B{ J)=8( J)-A { J,I) *T{I) 

30  CONTINUE 

00  35  1  =  1, N 
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35 


37 

41 

46 

50 

51 

52 
54 
62 
67 
67 


K=N+1-I 

B(K)=TCK)/A (K,K) 
K1=K-1 

IF(K1.£Q.0)  GO  TO  35 
DO  35  J1=1,K1 
J=K-J1 

T(J)=T( J)-A( J,K) *B(K) 

CONTINUE 

RETURN 

END 


6 

6 

7 

10 

11 

12 

13 

20 

21 

23 

25 

32 

33 

36 

37 
41 
47 
47 


FUNCTION  FU(I,A,B) 

COMMON/AUX/H,P,PKl,PK2,EMU,X,Y 

X=A 

V=8 

IF{A*B)5,10,5 
10  FU=0.0 

RETURN 

5  SUH=SIMF(I,fl.0,5.0) 

£R=0.01 
DEL  =5.  0 
20  UP=D£L+5.0 

AOOL=SIMF(I,D£L,UP) 

DEL  =UP 

T£ST=ABS(AODL/SUM) 

SUN=SUM+ADDL 
IF  (TEST-ER) 15,20,20 
15  FU=SQRT {X*Y) ♦SUM 
RETURN 
END 


3 

3 

5 

6 
15 
24 
31 
31 

33 

34 

36 

37 

41 

42 

44 

45 
62 

62 

63 


SU6ROUTINE  CONST (I) 

COMMON/ AUX/H,P,PK1,PK2,BMU,X,Y 

PR1=0,2S 

PR  2  =  0,29 

PK1=SQRT( (1. -2. »FR1) /{ 2.^(1. -PR  1)  )  ) 

PK2  =  SQRT(  (l.-2.*PR2)  /  (2 (1 , -PR2)  )  > 

READ  1,P 

1  FORMATCF10.5) 

HH=0.1 
HH=1Q. 0 
HH=5. 0 

HH=4. 0 
HH=1.0 
HH=0, 5 
HH=2. 0 
H=1,/HH 

PRINT  2,BMU,PR1,PR2,HH,F 

2  FORMAT(/////5X,*  MU2/MU1  =*F6.2,*  NUl  =»F4.2,»  NU2  =^F4 . 2/// 5X A 
1/H  =»F4.2,*  C21/PA  =*Fk.2) 

RETURN 

END 


5 

5 

23 

25 

27 

30 

31 
31 

33 

34 
36 
!fO 
L6 
55 
64 
73 
77 


FUNCTION  Z(I,S) 

COMMON/ AUX/H ,P, P K1 ,P K2 , EMU , X , Y 
B£SJH(A)=SQRT(2.*A/FI) *SIN( A) /A 
PI=3. 1415926 
IF(S-0,  0)5,5,10 
5  Z=0.0 
RETURN 
10  CONTINUE 

pp_p»p 

C1=PK1»PK1 
C2=PK2*PK2 
CC= 1 , ”C 1 

GA=SQRT(S»S+C1/PP) 

G3=SQRTCS»S+1./PP) 

GC=SQRT(S*S+C2/BMU/FF) 

GD=SQRT (S*S+1,/BMU/PP) 
AA=S»S+l./PP/2. 

A3=l. -3MI 
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100 
103 
110 
115 
123 
131 
137 
145 
152 
157 
162 
164 
171 
174 
2  00 
203 
206 
211 
214 
217 
222 
226 
231 
233 

235 

236 
240 
242 
244 

246 

247 
252 
256 
262 
266 
275 
301 
307 
312 
321 
332 
337 
347 
35  0 


Sg;fGi:iEr/Sc/PP/2.;BMu 

A6=  CS*S+GA»GOI  /AC/PP/2.  *EMU 
AHsfS*S-GB*GC)/AC/PP/2.*EHU 
AI=(S*S+GB*GC>/AC/PP/2.*BMU 
AJ=(GA-GC)/AC/PP/2.’0MU 

AK^CGA+GC) /AC/PP/2»*BMU 
A1=-(AB*GB-A0I 
A2=AB*GB-A|^.  _ 

A3sAA-BMU*S*S-AF 
A4= AA”fiMU^S*S~AG 
A5=-AA;eMU*S*S+AH 
A6=-AA*BMU»S*S+AI 
A7=S*CAE*GA-AJ) 

A8=-S»(AE*GA-AK) 

BA=A1*A6-A2*A5 
BB=A3*A6-S*A2*A7 
BC=A4*A6-S*A2*A8 
BD=S*A1^A7-A3*A5 
BE*S^A1*A8-A4*A5 
-ai»8B/BA - 

B2=BC/8A 
B3S8D/BA 
04seE/BA 
EA=2.*6A*H  - 

EBs2.*G B*H 
lc={EA+EE)/2. 

ED=2.*EC^^^ 

ElseXP(-EA)  - 

DL=B2>B3*E44'B4*£2+B1*E1 

Dl=2.*PP/CC/GE/OL_ 

D2=AA*AA;S*S*GA*GB 

DSI?*AA*(68»{Bl*e4-B2*e3>-S*S*GA>*E3 

D5=(AA*AA+S*S»GA*GB)*(B4*E2-91*£1) 

|;?F‘-l?S'eii3Sl!sS5’)-BESJH(S‘Y) 

RETURN 

END 


inversion  integral 

OIHENSION  »C<4,50) 

DIMENSION  3K(101)»TTC101) 

COMMON/ 2/TI*TF,DT ,MN,EK»TT 
READ  1,NN,MN,MM 

1  FORMAT(3I2) 

READ  2»7IjTF»pT 

2  FORMAT(3F10*5) 

PRINT  9S 

99  FORMATtlHl)  _ut  vm  rs 

GALL  SPLICE CGLAM^PHI fKMfC) 

102  FORMAT(5X,F10.5>5X,F10. 5) 

M11*MM-1 

300  F0RMAT(////5X,*  CC1*J)  C(2«J)  CC 

SS^«  121,  HN 

READ  3 ♦BET ,D£L 

3  FORMAT (2F1 0.5) 

PRINT  9  8*BET,Dc.L 


DELTA  =+F5.3) 


140 

140 

143 

144 
150 
153 
155 
161 
165 
165 
175 
175 
175 
177 
201 
202 

203 

204 
206 
210 
212 

214 

215 

216 
221 
223 

225 

226 
227 
230 
233 
235 
235 
237 
237 
237 
241 

244 

245 

246 
250 
252 
252 
254 
256 
26  0 
261 
264 
266 
27  0 
273 
275 
301 
304 

306 

307 
313 
315 

320 

321 
325 

325 

326 


98  F0RHAT(/////5X,*BETA  =»F5.3,» 

DO  11  L=1,MN 
AL  =  L 

S=1./<AL»BET)/D£L  ^ 

CALL  SPLIN£CGLAM,PHI ,MM,C»S,G) 

F=G»S 

IF(AL“2.)81, 82,83  ^ 

81  AC1)=(1.+BET)^0£L»F 

82  A(  2l=(H.*BET)*0EL»F-A(l))*  (3.  + BET) 

GO  TO  11  - 

83  CONTINUE 
T0P=1. 

L1=L-1 
AL1=L1  - 

DO  12  J=1,L1 
AJ=J 

TOP=AJ+TOP 

12  CONTINUE 
L2=2*L-1 
aoT=i. 

DO  13  J=L,L2 
AJ=J 

E0T=(AJ+EET)’B0T 

13  CONTINUE 
MUL=eOT/TOP 
SUM=0.Q 

DO  14  N=1,L1 
AN=N 

IF(AN-2.)  85,86,87 

85  TOC-1. 

GO  TO  88 

86  TOD=ALl 
GO  TO  88 

87  CONTINUE 
TOD=l. 

ICH  =  L1-  tN-2) 

DO  15  J=IGH,L1 
AJ=J 

TOC=AJ»TCO 

15  CONTINUE 

88  CONTINUE 
80D=1. 

JA=L1+N 

DO  16  J=L,JA 
AJ=J 

aOD=eOO^ (AJ+BET) 

16  CONTINUE 
CO=TOD/BOD 
SUM=SUN+CO  +  A  (N) 

14  CONTINUE  .  ^ 
A(L)=MUL»{OEL^F-SUM) 

11  CONTINUE 

CALL  JACSERCOEL,A,BET) 

CALL  NAMPLT  .  „  „  „  „ 

CALL  QIKSET(6.0, O.D, 0.0,e.0,0.0,0. 
call  QIKSAX(3,3> 

CALL  QIKFLT  (TT,eK,101) 

CALL  ENDFLT 
10  CONTINUE 
999  CONTINUE 
RETURN 
END 


6 

6 

6 

6 

7 

10 

11 

12 

14 

24 
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SUBROUTINE  JAGSER (D, C ,8)  ^  ^ 

DIMENSICN  CC50)  ,SF4|0)  ,P  (50) 

DIMENSION  BKdOl)  ,TT  (101) 

COMMON/ 2/TI ,7F,DT,MN,BK,TT 

TT(1)=0.0 

BK(1)=0.0 

LH=1 

T=TI 

T=T+OT 

X=2.*EXP (-D*T)-1. 

GALL  JACOBI (MN,X,e,P) 
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0) 


SF(1)=C(1)^P(1) 

DO  10  L=2,MN 

L1=L-1 

ALsL 

SF(L)=SF(L1) +CCL1*P(L) 

10  CONTINUE 

PRINT  97,T,X  .  ^ 

97  FORMAT(/////5X,»  T  =*F6.3,*  X 

DOT  |JT  Qg 

96  FORMAT(///5X,*  I  CCI)  *»5X, 

00  11  1=1.6 

PRINT  95»I,C(I>,I»SFCI) 

95  FORMAT (5X, 12, FI 0.2. 5X,I2.F1 0.5) 

11  CONTINUE 
LM=LM+1 
BK(LM)=SF(5) 

TT(LM)=T 

IF(T.LE.TF)  GO  TO  12 

RETURN 

END 


=-*F10.5) 

‘  N  F(T) 


THIS  PRCgIaM  CALCUtiT^S  JACOBI  POLYNOMIALS  OF  CROtR 
K-1  WITH  ARG  X  AND  PARAMETER  B  GT  -1 
OIHENSIC)  PB(N) 

AN>N 

IF(AN-2.)1,2,3 

P3(l)=l. 

RETURN  . 

Pe{l)=l. 

PB(2)=X-e*  (l.-X)/2. 

RETURN 

aSQ=8*B 

80NE=B+1. 

PB(1)=1. 

P3{2)=X-S*(l.-X)/2. 

DO  A  K=3,N 
AK=K 

AK1=AK-1. 

AK2=AK-2. 

K1=K-1 

K2=K-2 

C01=({2.*AK1)+8)*X 
C01=C (2.»AK2)+9)*C01 
C01=C<2.»AK2) ♦90NE)*(C01-8SQ) 

C02=2.*AK2*CAK2  +  3) (2.»AK1)  +B) 

C0=2.*AK1*  (AKH-B)»({  2.»AK2)  +B) 
PB(K)=(C01*P6(K1)-C02*PB (K2) )/C0 
RETURN 
END 


SUBROUTINE  SPLINE <X , Y .M ,C , XINT ,YINT ) 

DIMENSION  X(50)  ♦  Y(50)  ,C(*i,50) 

IFJXINT-XD)  1,10,11 
YINT=Y{1) 

RETURN 

CONTINUE 

IF(X(M)-XINT)1,12,13 

YINT=V(M) 

RETURN 

CONTINUE 

K=M/2 

N=M 

CONTINUE 

IF(X(K)-XlNT)3,li|,5 

YINT=Y(K) 

RETURN 

CONTINUE 

IF(XINT-X  (K+D)  A,  15,7 
YINT=Y(K+1) 

RETURN 

YINT=  (X  (K  +  l)  -XINT)  ♦(C(1,K)  +  (X(K+1)  -  XINT  )  •*♦  2  +  C  ( 3  ,K)  ) 
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54 

65 

65 

65 

70 

72 

72 

74 

75 
77 

100 

100 

102 
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106 

106 

111 

113 


114 

121 


c 


121 


123 


YINT=YINT  +  (XINT-X(K)  )  ♦  ( C  (  2  ,  K)  ♦  (  X INT-X  ( K )  )  ♦♦  2+ C  ( 4 ,  K)  ) 
RETURN 

5  CONTINUE 

IF(X(K-1>-XINT)6»16,17 

6  K=K-1 
GO  TO  4 

16  YINT=Y(K-1) 

RETURN 

17  N=K 
K=K/2 

GO  TO  2 -  : 

7  LL  =  K 
K=(N+K)  /2 

8  CONTINUE 

IF(X(K) -XINT)3,14,18 

18  CONTINUE 

IFCX(K-l)-XINT)  6»ie»19 

19  N=K 
K=(LL+K)/2 
GO  TO  8 

1  PRINT  101 

101  FORMATt*  OUT  OF  RANGE  FOR  INTERPOLATION  *) 

STOP 

END 


SUaROUTINE  SPLICECX,Y,M,C) 

DIMENSION  X(50)  ,Y(50)  ,0(5  0)  ,P(5fl)  ,E(50)  ,C(4,50) 
OIMENSICN  A (50, 3) ,a( 50)  ,Z(50) 

MM=M-1 
00  2  K=1,MM 
D<K)=X(K*1) -X(K) 

P(K)=D(K)/6, 

2  ECK)=(Y (K+13-Y(K)) /D (K) 

DO  3  K=2,MM 

3  B(K)=£( K)-E(K-l) 

A(  1,2)  =  -1.-D(1)/D(  2) 

A (1,3)  =  D (1) /D  (2) 

A{2,3)=P(2)-P(1)*A{1,3) 

A(2,2)  =  2.*(P(l)+°(2))-P(l)^A{l,2) 

A(2,3)  =  A(2,3)/A(2,2) 
a(2)  =  3(2)/A(2,2) 

DO  4  K=3,MH 

A(K,2)  =  2.*(P(K-1)  +P(K))  -F  (K-1)  ♦ACK-l,  3) 

B(K)  =  8(K)-P  (K-1)*B{K-1) 

A(K,3)=P(K)/A(K,2) 

4  a(K)  =  B(K)/A(K,2) 

0=0  (M-2)/D  (M-1) 

A{M,l)  =  l.+Q  +  A(M-2,3) 

A(M,2)=-Q-A (M,1)»A(M-1,3) 
e(M)=B(M-2> -A(M,l)»B(M-l) 

Z(M)=BCM)/A (H,2) 

MN=M-2 
DO  6  1=1, MN 
K  —  M  *^1 

6  Z(K)=9(K)-A(K,3)»Z  (Ktl) 

Z(  1)=-A(1,2)*Z(2)-A(1,3)*Z(3) 

DO  7  K=1,MM 
Q=1./(6.^0(K)) 

C(1,K)  =  Z(K)-*Q 

C{2,K)=Z(K+1)*Q 

C(3,K)=Y(K)/D(K)-Z(K)*PCK> 

7  C(4,K)=Y (K+1)/0(K)-Z(K+1)^P(K) 

RETURN 

END 


TO'Uyioncit  'iinpacZ 

PROGRAM  a£TA(INPUT,CUTPUT»PUNCH, PLOT, TAPE  99=PL0TJ 
REAL  NON (4) ,F(4,4 , 1) ,GC4,4)  ,  0(4) ,PT(4) 

REAL  B(4) ,C{4) 

REAL  LP(50) ,OTA C50) 

EQUIVALENCE  (NON,  3) 

_  COMMON  K1,K2,K3,K4  _ 

COMMON/A UX/H,P,PK1,PK2,BMU,X,Y 
LP(l)-0. 0 
DTA(1)=0 .0 
READ  2,K1,K2,K3,K4 
2  F0RMAT(I2) 

♦  <1  =  ORDER  OF  SYSTEM  OF  EQUATIONS 
»  K2  =  NO,  OF  DISTINCT  KERNELS 

♦  K3  =  NO.  OF  DATA  POINTS 

♦  K4  =  NO.  OF  DATA  SETS  TO  BE  EVALUATED 

♦  SET  UP  DATA  POINTS 

AK=K3 

DO  5  N=1,K3 
A  tcN 

5  PT(N)=AN/AK 

♦  SET  UP  INTEGRATION  MATRIX 

M=K3-2 

N=:K3-1 

A=K3 

A=1./C3.*A) 

DO  1(5  K=2,M,2 
10  0(KI=2.»A 

DO  15  K=1,N,2 
15  D(K)=4.*A 

0  (K3)sA 

♦  CALCULATE  NONHOMOGEN  ECUS  TERMS 

RHSsl.O 
00  22  1=1, K2 
PRINT  9 
9  FORMATdHl) 

READ  61,BMU 
61  FORMATfFlO. 5) 

DO  999  11=1, K4 
D  C  35  N=1,K3 

35  NON(N)=RHS*PT(N)*PT(N) 

♦  CALCULATE  KERNEL  MATRICES 

CALL  CONST (!» 

DO  20  N=1,K3 
DO  20  M=1,K3 
IF(M-N)25,30,30 
25  F  (M,N,I)=F(N,M,I) 

G  C  TO  20 

30  F CM,N,I) =FU(I,PT{H),PT(N») 

20  CONTINUE 

CALL  CHANGE(F,G,D,I) 

CALL  LINEQ(6,3,C,  K3I 
DO  40  L=1,K3 
PRINT  6,PT(L) ,NON (L) 

6  FORMAT(5X,F8.4,F15.6) 

40  CONTINUE 

LPdH-l)  =N0N(K3I 
OTA(II-H)=P 
QQQ  CONTINUE 

PUNCH  66,C0TA(IX»  ,  LP  ( IX>  ,IX  =  1,  19  ) 

66  FORMAT(2F10.5) 

CALL  LAPINV(OTA,LP) 

22  CONTINUE 
END 


FUNCTION  SIMP(I,A,B) 

COMM  ON/A  UX/H,P,PK1,PK2,BMU,X,Y 
OEL=0.25*(3-A) 

IF(OEL)40,45,50 

SIMP=0.0 

RETURN 

CONTINUE 

S A=Z  (I,A )+Z  (1,31 
SB=Z(I,A+2.»OEL) 

SC=Z  (I,A  +  0ELI+Z(I  ,A+3.*DEL) 


Sl=(0EL/3.)  *(SA+2.»SB+4.»SC) 

IFCSl.EQ.O. 0)  GO  TO  45 
K=d 

35  S8=SB<-SC 

D  EL  =  0,5*DEL 
SC=Z  {I*A  +  DEL) 

J  =  K-1 

00  5  N=3,J,2 
AN=N 

5  S  C=SC»-Z(I,A+AN»OEL) 

S2=  ( DEL/  3  . )  *  ( SA  +2  .  ♦S8+4,  *SC ) 

DIF=ABSl (S2-S1)/S1> 

ER=0.01 

IFIDIF-ER)30,25,25 
30  SIMP=S2 
RETURN 
25  K=2*K 

S1=S2 

IF(K-2043)35,35,40 
40  PRINT  42, I, A, 8 

42  FCRMAT(5X,»  INT.  DOES  NOT  CONVERGE  *,I3,2F9.4) 
PRINT  60,X,Y 
60  FORMAT(2F10.5) 

DO  70  J=l,10 
DIP=J 

D IP=0IP/lfl. 

W=Z(I,DIP) 

PRINT  60, W 
70  CONTINUE 
CALL  EXIT 
END 


SUBROUTINE  CH ANGE  (F ,G,D,  I) 
7  REAL  F<4,4, 1) ,G (4 ,4) ,D (4) 

7  COMMON  K1,K2,K3,K4 

7  DO  10  N=1,K3 

1C  DC  10  M=1,K3 

11  G(M,N)  =F(M,N,I)  ♦D(N) 

22*  10  CONTINUE 

30  DO  20  N=1,K3 

31  20  G  »N,N)=G(N,N)<-1.0 

40  RETURN 

41  E  NO 


7 

SUBROUTINE  LINEQ(A,3, 
REAL  A(N,N)  ,3(N),  T(N) 

7 

DO  5  1=2, N 

10 

5 

A  (I,1)  =  A  (I, D/A  (1,1) 

17 

DO  10  K=2,N 

20 

M=K-1 

22 

DO  15  1=1, N 

23 

15 

T(I)=A(I,K) 

33 

DC  2  0  J=1,M 

34 

A (J,K)=T (J) 

41 

J1=J+1 

43 

0  C  20  I=J1,N 

44 

T(I)=T(I)-A(I,J)*A(J, 

55 

20 

CONTINUE 

€1 

A (K,K)=T (K) 

€5 

IF(K.EQ.N)  GO  TO  10 

(6 

M=K<-1 

10 

00  25  I=M,N 

'll 

25 

A  (I.K)=T(I)/A(K,K) 

1  '5 

10 

CONTINUE 

»  BACK  SUBSTITUTE 

1  jO 

DO  30  1=1, N 

1  il 

T (I)=B(I) 

1  14 

M  =  I  +  1 

1  re 

IF(M.GT.N)  GO  TO  30 

1  -'1 

DO  30  J=M,N 

1  12 

B  {J)=8(J)-A(J,I)*T(I) 

1  '2 

30 

CONTINUE 

1  ^6 

DO  35  1=1, N 
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.  Z7 
L  tl 
.46 
.50 
LSI 
.52 
L54 
.62 
L67 
.67 


3lK)=t^K)/A (K.K) 

1/  ^  s  1C**  i 

IF(Kl.EO.fl)  GO  TO  35 
DO  35  J1=1*K1 

T  uT=TIJ)-A(J,K)*BCK) 
35  CONTINUE 
RETURN 
END 


6 

6 

7 

10 

11 

12 

13 

» 

21 

25 

c2 

:3 

‘6 

Z7 

n 

n 

n 


10 

5 

20 


FUNCTION 

COMMON/A UX/H,P*PKl,PK2*aMU*X,T 

XsA 


Y 

IF(A»a)5  tl0»5 

FU=0.0 

RETURN 

SUMsSIMP  C I»0*0*  5. 
ERs0*01 
DEL  =5.0 
UP=0EL+5 • 0 
AOOL=hMP(I,DEL,U 
D  EL  =UP 

TEST=ABS (ADDL/SUH 
SUH=SUH+AOOL 
IF(TEST-ER) 15,20, 
FU=SQRT<X*Y)*SUH 

RETURN 

END 


0) 

P) 

) 

20 


3 

3 

5 

6 
15 
24 
31 
31 

33 

34 

36 

37 

41 

42 
57 

57 

60 


PRlffl.29 

oeF-Eni?  ((l.-2.*PRll  /  (2.*  (l.-PRl  >  J  * 

PK2=iQRT { (i.-2.*PR2)/C2.»(l.-PR2>)> 

READ  1,P^^ 

FORMAT(F10.5I 
H  Fs5.0 
HH®0.2 
HH=0.5 
HH=1.0 

NUl  =*F4.2,*  NU2 

1/H  =*F4.2,*  C21/PA  =*F4.2) 

return 

END 


==F4.2///5X,*  A 


5 

5 

26 

30 

32 

33 

34 
34 
36 
45 

-  54 
60 
62 
72 
.02 
05 
.16 
17 


Z=0.D 
R  ETURN 
CONTINUE 
ppsp*p 

GB=SORTtS*S  +  l./PP) 

AO=l •♦AA/AB^EXP C -2.*GB*H) 
Z=^F-sf*3ES  JT  {S*X  )  *BESJT  (S*Y) 

return 

END 
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5 

5 

5 

5 

5 

16 
16 
30 
30 
34 
34 
40 
44 
44 
65 
65 
67 
7  3 

73 

12 

12 

j6 

?1 

30 

30 

40 

40 

43 

44 
50 
53 
55 

a 

f5 

•5 

5 

'5 

75 
77 
fl 
i2 
^3 
'4 
C6 
jO 
i2 
j4 
!5 

:s 

71 

'3 

25 

76 
27 

:o 

23 

25 

‘5 

27 

27 

27 

il 

t4 

5 

i6 

20 

-22 

to 

"4 

‘:6 

«0 


SUBROUTINE  LAPI NV ( GLAK,PHI ) 

C  THIS  PROGRAM  EVALUATES  THE  COEFFICIENTS  FOR  SERIES 

C  OF  JACOBI  POLYNOMIALS  WHICH  REPRESENTS  A  LAPLACE 

C  INVERSION  INTEGRAL 

REAL  MUL 

DIMENSION  A  (50)  ,G LAM (50)  , PHI  (5  0)  ,C  (4,50) 

DIMENSION  BKClOl) ,TT(101) 

COMMON/2/TI,TF,OT ,MN,BK,TT 
READ  1,NN,MN,MH 

1  F0RMAT(3I2) 

READ  2,TI,TF,0T 

2  FORMAT(3F10.5) 

PRINT  99 

99  FORMAT(lHl) 

CALL  SPL ICE (GLAM,  PHI, MM, G) 

PRINT  101 

101  F0RMAT(/////5X,*  GLAM  PHI  *) 

PRINT  102,  (GLAMd  )  ,PHI(I),I  =  1,MM) 

102  FORMAT(5X,F10.5,5X,F10.5  ) 

M11=MM-1 

PRINT  300 

303  F0RMAT(////5X,*  C{1,J)  C(2,J)  C(3,J) 

1,J)  *) 

PRINT  103,  (  {C(I,J)  ,I=1,4),J  =  1,  Mil) 

10  3  F  CRHAT(5X,F1C  .5,5X,F10,5  ,5X,F10.5,5X,F10.5) 

PRINT  99 
DO  10  1=1, NN 
READ  3, BET, DEL 

3  FORMAT(2FlQ.5) 

PRINT  98, BET, DEL 

93  F0RMAT(/////5X,*aETA  =*F5.3,*  DELTA  =»F5.3) 

DO  11  L=1,MN 
A  L=L 

S=l. /(AL+BET) /DEL 

CALL  SPLINE  (GLAM,  PHI, MM, G,S,G) 

F=G*S 

I F(AL-2. )81,82, 83 
81  A  (1)  =  (1.*BET)  »DEL*F 
GO  TO  11 

S2  A (2)=((2.+BET)*0EL*F-A{l))*{3,+BET) 

GO  TO  11 
83  CONTINUE 
TOP=l. 

L1=L-1 
AL1  =  L1 

DO  12  J=1,L1 
A  .=  J 

TOP=AJ*TOP 

12  CONTINUE 
L2=2»L-1 
80T=1, 

DO  13  J=L,L2 
A  J=J 

3OT=(AJ+3ET)*30T 

13  CONTINUE 
MUL=80T/T0P 
S  LM  =  0.G 

DO  14  N=1,L1 
A  N=N 

IF(AN-2.  )  85,86,  87 
85  TOO=l. 

GO  TO  88 
36  TOD=ALl 
GO  TO  33 
87  CONTINUE 
TOD=l. 

ICH=Ll-(N-2) 

DO  15  J=ICH,L1 
A  J=J 

TOD=AJ*TOD 
15  CONTINUE 
33  CONTINUE 
BOO=l. 
jA=LH-N 
DO  16  J=L,JA 
A  J=J 


C{4 
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BOO=BOO* (AJ+QST) 

CONTINUE 

CO=TOD/BOD 

SUM=SUM+CO*AiN» 

CONTINUE 

A (L)=MUL*(DEL*F-SUN) 

CONTINUE-- 

C  ALL  JACSER(D£L,A  ,BET) 

CALL  NAMPLT  ^ 

CALL  QIKSET(6.0.0.  C,0.0, 6,  0,  0.  0*  0. 0> 
-CALL.QIKSAX{3t3) 

CALL  QIKPLT(TT,3K,101) 

CALL  ENDPLT 

CONTINUE 

CONTINUE 

RETURN 

END 


SUBROUTINE  JACSER  (0,0,8) 

DIMENSION  C(50)  ,SF(50),P(50) 

DIMENSION  BKIlOl)  ,TT{101) 

COMMON/2 /TI,TF,0T,MN,8K,TT 

TT(1)=0,0 

3K(1)=0. 0 

LM=1 

T=TI 

T=T+OT 

Xs2.*EXP(-0*T)-l. 

CALL  JAC03ICMN, X,  3,P) 
SFtl)=C(l)*P(l) 

DO  10  L=2,MN 

L1=L-1 

AL=L 

SFCL)=SF(L1) tCtL) *P(L) 

CONTINUE 
PRINT  97,T,X 

FORMAT(/////5X,  ♦  T  =*F6.3,*  .X 
PRINT  96 

F0RMAT(///5X,*  I  C(I)  ♦♦SX, 
DO  11  1=1,6 

PRINT  95,1, GCI), I, SF(I) 

FORMAT(5X,I2,F10. 2,5X,I2,F10 .5 ) 

CONTINUE 

LM=LM+1 

BK(LM)=SF(5) 

TT(LM)=T 

IFfT.LE.TF)  GO  TO  12 

RETURN 

END 


=*F10.5} 

N  F(T) 


SUBROUTINE  JACOBI  ( N,X,3,P3) 

THIS  PROGRAM  CALCULATES  JACOBI  POLYNOMIALS  OF  ORDER 
K-1  WITH  ARG  X  AND  PARAMETER  8  GT  -1 
0  IMENSION  P8(N) 

A  )=N 

I  F(AN-2.  )1,2,3 
PB(1)=1« 

RETURN 

P3(l)=l. 

Pa(2)=X-B*(l.-X)/2. 

RETURN 

BSQ=B*B 

30NE=3<-1. 

PB<1)=1. 

P8(2)=X-B»{1.-X)/ 2. 

DO  4  K=3  ,N 
AK=K 

AK1=AK-1. 

A  K2  =  AK-2. 

K1=K-1 

K2=K-2 

C01=((2.*AK1)«-9)*X 
C01=  ( (2,  ♦AKE)  4-8)  ’COl 


51 

C01= ( (2. *AK2) +80NE  )*CC01-BSQ) 

56 

C02  =  2.»AK2*  CAK2  +  B )♦< {2.»AK1>  fB) 

64 

C0=2.*AK1*(AK1+B)  ♦  (  {2.»AK2) +B) 

71 

4 

PB(K)  =  (C01*PB(K1)  -C02*PB(K2) )/CO 

1  02 

RETURN 

103 

END 

SUBROUTINE  SPLINE  (X,Y,M,C,  XINT,YINT) 

11 

DIMENSION  X(501  tY  (50),C(4,50) 

11 

IF(XINT-X(1))1,10,11 

13 

10 

YINT=Y{1) 

14 

RETURN 

15 

11 

CONTINUE 

t5 

IF  CX  (M)-XINT)  1,  12  »  13 

;i 

12 

YINT=Y(M) 

-3 

RETURN 

'.3 

13 

CONTINUE 

?3 

K=M/2 

^5 

N=t1 

76 

2 

CCNTINUE 

76 

IF{X{K)-XINT)  3,14 ,5 

;2 

I'* 

YINT=Y(K) 

4 

RETURN 

75 

3 

CONTINUE 

:5 

IF{XINT-X(K  +  1})4,  15,7 

J-l 

15 

YINT=Y(K+1) 

13 

RETURN 

’-3 

4 

CONTINUE 

-3 

YINT={XCK  +  1)-XINT  )*{CC1,K)»(X(K+1)-XINT)  **2  +  C(3, 

54 

YINT=YINT+(XINT-X  (  K )  )  ♦  (C  (2  , K  )  ♦  (X  INT- X ( K )  >  ♦♦  2 +C  ( 4 

65 

RETURN 

65 

5 

CONTINUE 

65 

IF(X  CK-1)-XINT)  6,  16,17 

70 

6 

K=K-1 

72 

GO  TO  4 

72 

16 

YINT=Y(K-1) 

7h 

RETURN 

75 

17 

N=K 

77 

K=K/2 

130 

GC  TO  2 

100 

7 

L  L=K 

132 

K=(N+K)/2 

103 

8 

CONTINUE 

103 

IF(X(K)-XINT) 3,14,13 

106 

18 

CONTINUE 

106 

IF(X<K-1)-XINT) 6,  16,19 

111 

19 

N=K 

113 

K  =  {LL+K)  /2 

114 

GC  TO  8 

1  15 

1 

PRINT  IGl 

1  71 

101 

FORMAT(»  OUT  OF  RANGE  FOR  INTERPOLATION  ♦) 

1  21 

STOP 

123 

END 

7 

SUBROUTINE  SPLI CE  (  X,  Y  ,M,  C) 

0  IMENSION  X  C50)  ,  Y  (50  )  ,0  (50)  ,P(50  )  ,E(50),  C(4,  50  ) 

7 

DIMENSION  A  (50,3)  ,8(50) ,Z(50) 

7 

MM=M-1 

11 

DO  2  K=1,MM 

12 

D(K)=X(K+1)-X(K) 

15 

P  {K)=0(K)/6. 

20 

2 

E  (K)  =  (Y(Ktl)-Y(K)  )  /D(K) 

26 

DO  3  K=2,MM 

27 

3 

B  (K)=ECK)-£(K-1) 

34 

A  (1,2)=-1.-0(1)  /D  (2) 

37 

A(1,3)=D (1)/D(2) 

41 

A  (2,3)=P  (2)-P(l)*  A  (1,3) 

44 

At2,2)=2,*(P(l)+P  (2)  )-P(l)*A(l,2) 

50 

A (2, 3)=A (2,3) /A (2 ,2) 

3  (2)=3(2) /A (2,2) 

51 

53 

DO  4  K=3,MM 

54 

A  (K,  2)=2.*  (P(K-l)  +P(K))-P(K-1)  *A  (K-l  ,3) 

61 

B  (K)=B(K)-P(K-i)*B(K-l) 
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A{Kt3)=P(K)/A(K,2l 
4  3  (K)=8(K )/A (K,2I 
Q=0(M-2) /0{M-1) 

A  {M,l)=l.+Q+A{M-2»3) 

A  (M,2)=-Q-A  (H,1)*A(M-1,3) 

3 (M)=B(M-2)-A(M,l )*B(M-1) 

Z  (M)=B(M)/A  (M,2) 

MN=M-2 
00  6  1=1 tMN 

6  Z(K)=8(K)-A  (K,3jl*Z(K+l)  . 

Z  tl)  =-»A(l,2)*Z(2)  -AClt3)  *2(3) 
00  7  K=1  *MM 

0=l,/(6.  *0{K)  ) 

C  {1*K)=Z  (K)  *Q 

C t2,KI=Z CK+1) *0 

C  (3tK)=Y  (K)/0(K»-Z(K)*P(K) 

7  C  (4,K)=Y  (K«>1)/D(K  )-Z(K«-l)*P(K) 
RETURN 

ENO 
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